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We investigate a recently developed approach^ that uses maximally-localized Wannier functions (MLWFs) 
to evaluate the van der Waals (vdW) contribution to the total energy of a system calculated with density- 
functional theory (DFT). We test it on a set of atomic and molecular dimers of increasing complexity (argon, 
methane, ethene, benzene, phthalocyanine, and copper phthalocyanine) and demonstrate that the method, 
as originally proposed, has a number of shortcomings that hamper its predictive power. In order to overcome 
these problems, we have developed and implemented a number of improvements to the method and show 
that these modifications give rise to calculated binding energies and equilibrium geometries that are in closer 
agreement to results of quantum-chemical coupled-cluster calculations. 
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I. INTRODUCTION 

Local and semi-local exchange-correlation functionals 
used in density-functional theory 3,4 (DFT) can not ac- 
count for the effect of long-ranged dispersion, or van der 
Waals (vdW), interactions. Dispersion interactions are 
crucial for weakly-bound systems, particularly where no 
covalent or ionic bonding is present, and often dominate 
intermolecular binding energies and equilibrium geome- 
tries. Incorporating vdW interactions in DFT remains a 
challenging task and a wide variety of methods have been 
developed, approaching the problem from many different 
perspectives^ — . In this work we focus on the method re- 
cently proposed by Silvestrelld^, which has been recently 

applied to various systemsi 4 . and implemented in a A " Maximally-Localized Wannier Functions 



in Sec. |TT]we recap the necessary background theory re- 
lating to MLWFs and Silvestrelli's method; in Sec. ITTT1 
we highlight some of the problems with the method as 
it stands, and describe our improvements; in Sec. ITVl we 
then present and discuss results for vdW-corrected total 
energies and equilibrium geometries obtained by applying 
these methods to a series of dimer systems and compare 
to quantum chemical coupled-cluster and semi-empirical 
vdW (DFT+D) approaches; finally, in Sec. [VI we draw 
our conclusions. 



II. THEORETICAL BACKGROUND 



number of modern electronic structure code o 18 i 19 . This 
approach uses maximally-localized Wannier functions^ 
(MLWFs) as a means of decomposing the electronic den- 
sity of the system into a set of localized but overlapping 
fragments, which may then be used to calculate a vdW 
correction to the DFT total energy by considering pair- 
wise interactions between density fragments as derived 
by Andersson, Langreth and Lundqvist^ (ALL). 

In this Article, we explore the parameters and approx- 
imations involved in Silvestrelli's method and improve its 
results where possible by modifying various aspects of the 
method. We apply the method and our proposed mod- 
ifications to a series of test systems, then to two more 
challenging systems, a phthalocyanine and a copper ph- 
thalocyanine dimer. We thus demonstrate that although 
this method can offer an easily implementable and com- 
putationally efficient way of calculating the dispersion 
correction to the energy with the possibility of improved 
accuracy (once some modifications are applied to it), it is 
largely dependent on a number of parameters and choices 
one can make. 

The remainder of the Article is organized as follows: 



Wannier functions^ are orthogonal localized functions 
that span the same space as the eigenstates of a single 
particle Hamiltonian. Consider the set of N occ occupied 
(valence) eigenstates {|^ m )} of a molecule. The total en- 
ergy is invariant with respect to unitary transformations 
among the eigenstates 



N occ 

n) = ^ ^ UjYin I Urn ) 



(i) 



If the unitary matrix U is chosen such that the result- 
ing A^occ orbitals {w n (r)} minimize their total quadratic 
spread, given by 

Q = J2 (K|r 2 K) - (w n \r\w n ) 2 ) = «r 2 >« " r 2 ) , 

n n 

(2) 

then they are said to be maximally-localized Wannier 
functions^ (MLWFs). Each MLWF is characterized by 
a value for its quadratic spread, 5^, and its centre, f n . 

In the construction of MLWFs it is sometimes useful to 
consider not only the valence manifold but also a range 
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of unoccupied eigenstates above the Fermi level — often 
those constituting the antibonding counterparts to the 
valence states. This not only allows the MLWFs to be 
more localize d 22 ! 23 but can also restore symmetries that 
would otherwise be broken arbitrarily through the con- 
struction of MLWFs for the valence manifold only. 

In order to do so, one defines an outer energy win- 
dow, £?win, consisting of N w [ n > N occ states, from which 
one may extract an optimal A^is-dimensional subspace 
(^win > ^Vdis > N occ ) using the disentanglement ap- 
proach described in Ref. 0, 



,opt\ 



E^; 

P =i 



dis I 



(3) 



where U dls is a rectangular N w { n x 7Vdi S unitary matrix. 
TVdis MLWFs may then be localized by suitable rotation 
of the optimal subspace in the usual manner: 



m=l 



or, in terms of the Bloch states: 



N dis AT win 
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(4) 



(5) 



Furthermore, an inner, or frozen, energy window may 
be defined if one wishes to make certain that a range 
of low-lying eigenstates are included in the optimal sub- 
space, for example, the occupied states. Algorithms for 
determining MLWFs from the eigenstates obtained from 
electronic structure calculations are implemented within 
the Wannier90 software package 25 . 

The single-particle density operator is given by 



A^occ 

E \u n )(u r , 



(6) 



It can also be written in terms of the N occ fully-occupied 
valence MLWFs, \w n ) or equivalent ly in terms of a larger 
set of TVdis disentangled MLWFs \w dls ) that span the oc- 
cupied subspace, which can be guaranteed by using a 
suitable frozen / inner window in the disentanglement pro- 
cedure, and that have occupancies 
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where we have substituted Eq. (pQ) and Eq. (|5]), respec- 
tively, into Eq. ©, and where the occupancies are given 
by 



N, 



fki — E 



c -^dis 
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We can write the density as a sum of diagonal (I = k) 
and off-diagonal (I ^ k) terms, 

A^dis N dis 

P(r) = E /« K S W| 2 + E /i>r dta (r)«C(r), 

1=1 l^m 

= Pd{v)+ P od{v). (10) 

It is important to note that in this form, p£>(r) alone 
integrates to the number of valence electrons 7V e , be- 
cause the mututal orthogonality of the MLWFs ensures 
|p D(r)dr = 0. 

In the case of considering MLWFs obtained from the 
manifold of occupied states only (iVdi S = N occ ), the the 
occupancy matrix is simply the identity matrix, fki = Ski > 
and the charge density in terms of the MLWFs is simply 
given by 



P(r) 



N occ 



Mr) | 2 



(11) 



It is worth noting that in the case of spin-degenerate 
systems, the occupancies must be scaled by a factor of 2. 

We have adapted the Wannier90 code to calculate the 
occupation matrices, and can choose to make a diago- 
nal approximation to the density by retaining only the 
first term of Eq. ([TQ|) . The effect of approximating the 
true density with the diagonal approximation will be dis- 
cussed later in Sec. IIVII in the context of the improve- 
ments, described in Sec. IIII| to Silvestrelli's method. 



B. Silvestrelli's method 

Silvestrelli's approach^ 2 - is based on the Andersson, 
Langreth and Lundqvist 7 (ALL) expression for the vdW 
energy in terms of pairwise interactions between density 
fragments p n (r) and pi(r f ), separated by a distance r n i, 



Evdw 
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where g n i(r n i) is a damping function 2 which screens the 
unphysical divergence of Eq. ([T2]) at short range, and 
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in atomic units. It should be noted that these expressions 
are only strictly valid in the limit of non-overlapping den- 
sity fragments. There are various forms for the damping 
functio n 26 ! 27 that might have a slight short-range effect 
but should not affect the long-range behaviour of the 
vdW energies. Here we chose to use the damping func- 
tion as proposed in the original paper by Silvestrelli 1 . 

Now, in accord with Eq. ([TT]) . the MLWFs obtained 
from the valence orbitals of a system provide a localized 
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decomposition of the electronic charge density, such that 
p n (v) = \w n (r)\ 2 , so that Eq. ([T3|) becomes 



32tT^ J| r |< rc J\r'\<r' c 



dr 



/ K(r)|K(r')| 



K(r)| + Mr')r 
(14) 

where r c is a suitably chosen cutoff radius obtained by 
equating the length scale for density change to the elec- 
tron gas screening length 2 ; we will revisit this point later. 

In order to make the calculation of the integrals more 
tractable, the charge density is approximated by replac- 
ing each MLWF w n (r) with a hydrogenic s-orbital that 
has the same centre r n and spread S n as the MLWF, and 
whose analytic form is is given by 



3 3/4 



\fltSn 
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-V3|r-f n |/S„ 



(15) 



which, on substitution into Eq. (fT4|) and after some alge- 
bra, gives 



o3/2 q3 

2 • 3 5 / 4 



F(S n , Si), 



where 



F(S n ,Si) 



dx 



dy 



x 2 y 2 e x e y 
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(16) 



(17) 



P = {Sn/Sif 2 , x c = VSr c /S n and y c = yfa'JSi. 
Whereas evaluating Eq. (jHj) using the true MLWFs re- 
quires a computationally demanding six-dimensional nu- 
merical integration, Eq. ([T7|) may be evaluated easily 
since it is only a two-dimensional integral that depends 
solely on the MLWF spreads and centres, not their de- 
tailed shapes or orientations. 

We note that in the case of a spin-degenerate system, 
since every MLWF is doubly occupied, the density of 
each fragment must be multiplied by a factor of 2 and, 
therefore, the C§ n \ integral in Eq. ([Tl]) must be scaled by 
a factor of y/2. 



III. IMPROVEMENTS TO SILVESTRELLI'S METHOD 



rotation of the valence subspace only, which arbitrarily 
break the symmetry (we will demonstrate examples of 
this phenomenon in Sec. IIV|). 



In order to account for the partial occupancy of the 
MLWFs, we make a slight modification to Silvestrelli's 
approach, explicitly introducing occupancies in the def- 
inition of the Cqui integral; since in the diagonal ap- 
proximation, the density of each fragment is now given 
by p n (v) = f™ n \w n (r)\ 2 , the expression for F(S n ,Si) in 
Eq. (p"T|) becomes 



F(S n , S t 



)= dx 

JO JO 



dy- 



x 2 y 2 e x e y 



c /(PVf%n)+e- v /Vfu 

(18) 

where the f™ n are given by Eq. (|9]). We will see in Sec.HVl 
that this seemingly simple idea can give rise to a marked 
improvement in the accuracy of the method. 



B. Modification to describe p-like states 

MLWFs describing only the valence manifold often take 
the form of well-localized functions centred on a bond be- 
tween two atoms, and are thus reasonably well-described 
by the approximation of replacing them with a suitable 
s-orbital. When anti-bonding states are included in the 
construction of the MLWFs, the resulting orbit als have 
more atomic-orbital character. This is demonstrated by 
the atom- centred p- like MLWF shown in Fig.JTJ It is clear 
that the density associated with such an MLWF will not 
be very well represented by a single s-like function at its 
centre. In order to approximate p-like orbitals appro- 
priately when calculating Cq, one could imagine using 
a suitably-oriented analytic expression for a hydrogenic 
p-orbital, for example, a canonical p^-orbital given by 



Pz(r) 



30 5 / 4 rcos6> 



/30r/2S 



(19) 



The approximations that go into the method described 
in the previous Section will clearly not always hold, and 
the need to examine them is clear. In this Section, we 
introduce our enhancements to the method that address 
possible drawbacks. 

A. Partly Occupied Wannier Functions 

Using a manifold of eigenstates that includes but is 
larger than the subspace spanned by just the valence 
states results in partly-occupied MLWFs that are gen- 
erally more localized and that better reflect the symme- 
tries of the system, as opposed to MLWFs obtained by 



which has been normalized such that its quadratic spread 
is (Pz\(y — r) 2 \p z ) — S 2 - As a consequence of the explicit 
angular dependence, using this function in Eq. ([T^]) would 
give rise to four-dimensional integrals for which analytic 
solutions are not readily available. Numerical evalua- 
tion of these integrals, for realistic systems, would be 
prohibitively computationally expensive. We solve this 
problem by identifying the p-like MLWFs in the system 
and replacing them with the hydrogenic form given in 
Eq. ([T9|) . Then, we further approximate each lobe (lower 
and upper) of this p-like orbital with two separate hydro- 
genic 5-orbitals of the form of Eq. ([T5]) . In order to do 
so, for each of the upper (+) and lower (— ) lobes of the 
orbital, it is necessary to know the spread S± and centre 
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Figure 1. Partly occupied p-like orbital on ethene molecule. 
In the method described here, each of the two lobes (coloured 
red and blue) is replaced by an s orbital and considered a 
separate fragment. 



r±, given by 

roo ^71-/2 r2ir 



(20) 



POO P7T / Z PZ7T 

S 2 ±= I / / r 4 p 2 (r) sin OdrdOdcj), 
Jo Jo Jo 

poo p2-k 

r± = r ± / / / r 3 cos#p 2 (r) sin OdrdOdcf) z, 
Jo Jo Jo 

(21) 



which, after some algebra, simplifies to 

75 



S± = 



8^2' 



, 155 A 
r± = r ± — = z, 



30 



(22) 
(23) 



where f and S are the original centre and spread, re- 
spectively, of the true MLWF. These expressions may be 
easily generalized to arbitrary orientations of the sym- 
metry axis of a p-like state by rotating the offset vectors 
(f ± — f ) accordingly. 

Thus, we have developed a formalism whereby the 
charge density due to MLWFs with p-like character can 
be represented by a pair of s-like hydrogenic orbitals with 
appropriate centres and spreads. In Sec. US we will show 
how this works in practice for calculating vdW energy 
corrections. 

In the relatively simple systems studied in this pa- 
per, the p-like orbitals are easily distinguished from other 
orbitals by their partial occupancies, given by Eq. ([9]), 
which are typically closer to 0.5 rather than 1. Alterna- 
tively, and especially for structurally more complex sys- 
tems, the shape of each MLWF could be characterized 
using the efficient method described in Appendix A of 
Ref. [28| as another means of automating the procedure of 
identifying p-like functions. 



C. Symmetry Considerations 

Minimizing the total spread Q with respect to the el- 
ements of the unitary matrix U, and thus producing 
MLWFs, has the effect of picking from the space of all 
possible unitary matrices one which produces the most 
localized Wannier functions accessible through optimiza- 
tion from a chosen initial guess. This is often enough to 



uniquely determine the MLWFs. In some cases, however, 
it does not give rise to a unique choice, even if the opti- 
mization procedure is perfect. For example, the atomic 
positions and electron density of the system may possess 
certain symmetry elements, such as rotations about a 
particular axis. Then there will exist a number of equally 
valid and degenerate representations of the MLWFs and 
their centres, which give the same spread, and are re- 
lated by symmetry. The minimization procedure breaks 
the symmetry by choosing one of these representations; in 
other words there will be a degree of arbitrariness in the 
final MLWFs. It is clear from Eq. ([T2]) that any degree of 
non- uniqueness of the centres will cause an undesirable 
variability of the vdW energy calculated in Silvestrelli's 
method. This is indeed what we observe in some of the 
examples below. Moving away from a description of the 
MLWFs using the valence states only, and towards using 
partly occupied MLWFs that include anti-bonding states 
and which retain the symmetries of the system, enables 
us to overcome these problems, as we demonstrate below. 



IV. APPLICATIONS 
A. Calculation Details 

For the application of Silvestrelli's method to the fol- 
lowing dimer systems we used the Quantum Espresso 
(QE) package 18 to perform the ground-state DFT cal- 
culations, and Wannier90 25 to obtain the centres and 
spreads of the MLWFs. Our results are compared to 
both the semi-empirical DFT+D metho d 29 ! 30 as imple- 
mented in QE, which is expected to give good asymptotic 
behaviour, and a wavefunction-based coupled-cluster 
approach, CCSD(T), which is considered the 'gold- 
standard' of quantum chemistry. 

The PBE 3 -^ generalized-gradient approximation for ex- 
change and correlation, except in the case of argon where 
the revPBE 3 - 2 - functional was used; norm-conserving 
pseudopotentials, and T-point sampling of the Brillouin 
zone were used throughout. We note that we have chosen 
to use revPBE for the argon system since PBE produces 
significant binding in rare gas dimers as it overestimates 
the long-range part of the exchange contributio n 12 ! 33 ! 34 . 
For all the other systems we studied in this manuscript, 
however, PBE does not cause spurious binding and would 
therefore normally be considered an appropriate func- 
tional. A plane-wave basis set cut-off energy of 80 Ry 
was used in all calculations with QE except for the case 
of the phthalocyanine and copper phthalocyanine where 
a 50 Ry energy cutoff was used. For the dimers of argon, 
methane, ethene, phthalocyanine and copper phthalocya- 
nine, cubic simulation cells of length 15.87 A, 15.87 A, 
21.16 A and 23.81 A, respectively, were used. For the 
dimers of benzene, a hexagonal cell with a = 15.87 A 
and c = 31.75 A was used. For all the systems, the 
choice of energy windows when using the disentangle- 
ment procedure in Wannier90 for our modified method 
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was as follows: inner (frozen) energy windows were cho- 
sen to include all the valence states; outer energy win- 
dows ranged from the lowest eigenvalue of the system, 
e , to a maximum of E win = e L uMO + ^(^homo - eo), 
where chomo is the energy of the highest occupied va- 
lence Kohn-Sham (KS) state and 6lumo is the energy 
of the lowest unoccupied KS state. The factor a = 0.4 
was chosen to scale down the valence energy bandwidth, 
used to estimate the energy difference required above the 
LUMO when including anti-bonding states. We discuss 
the sensitivity of the method to this factor in Sec. IIV Jl 



B. Argon 

We will first investigate the severity of the aforemen- 
tioned issues relating to symmetry, by considering the 
case of an argon dimer. Optimization of the MLWFs 
describing a single argon atom produces four doubly oc- 
cupied MLWFs arranged tetrahedrally around the atom. 
Due to spherical symmetry, the orientation of these ML- 
WFs with respect to a given coordinate system is arbi- 
trary for an isolated atom and the final MLWFs obtained 
will depend on the initial guess used. In the dimer, this 
arbitrariness is removed, at least in principle, since the 
spherical symmetry is broken by the presence of the other 
atom at a specific orientation. At large separations, this 
is not in practice necessarily the case: the electron den- 
sity overlap between the Ar atoms is vanishingly small, 
since the wavefunctions decay exponentially away from 
the atom. Therefore, to within attainable numerical pre- 
cision, the orientation of the MLWFs on each atom is 
uncorrelated with the orientation of the other atom: the 
MLWFs can be freely rotated with respect to the atom 
without affecting the total spread. Note, however, that 
since the vdW energy only decays as R~ 6 , its value is 
influenced by the orientation of the MLWF centres (and 
hence their separation) out to distances beyond which 
the calculated spread (and thus the optimised MLWF 
orientation) has ceased to be sensitive to separation. 

This dependence can be investigated in a two-atom sys- 
tem by fixing the relative orientations of the MLWF cen- 
tres between the two atoms in the dimer. This is achieved 
by first calculating the MLWF centres for a single atom 
of argon and then translating and rotating these centres 
to the second Ar atom with various choices of alignment. 
We will refer to this approach as the fragment method. 
In this method, we calculate the dispersion correction to 
the energy for a dimer system using various possible ar- 
rangements of MLWF centres on the other atom. Three 
possible high-symmetry choices are shown in Fig. [2j For 
each of these orientations, Fig. [3] (top) shows the binding 
energy of the Ar dimer as the separation of the atoms 
varies. We see that there is considerable displacement of 
the curves, and the binding energy and the equilibrium 
separation change according to the alignment chosen by 
up to 0.04 kcal/mol and 0.08 A, respectively. 

In contrast to this fragment approach, in Fig. [3] (bot- 
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Figure 2. Illustration of three of the many possible configura- 
tions of MLWF centres (small pink spheres) for the two argon 
atoms (large blue spheres) in the fragment method. 



torn) we show the binding energy as calculated with the 
normal approach of using the optimized MLWFs of the 
entire dimer system. However, here we have used vary- 
ing initial guesses corresponding to the set of possible 
alignments shown in Fig. [2j We see that at small sepa- 
rations, the MLWF centres always converge to the same 
positions, regardless of the initial guess, and the binding 
energy curve is nearly independent of the choice of initial 
guess (~ 10 ~ 3 kcal/mol variation). 

At larger separation, however, the spread minimization 
is insufficiently sensitive to the relative orientation of the 
MLWFs on different atoms, and does not necessarily al- 
ter it from the initial guess, resulting in several different 
possible results depending on the initial orientation of 
the centres. If a random initial guess is chosen, then 
the energy varies discontinuously, as a function of sepa- 
ration, within the bounds imposed by the limiting cases 
described using the fragment method. This is because the 
MLWF centres converge to different orientations depend- 
ing on their starting positions (curve labelled 'random' 
in Fig. [3] (bottom)). 

In order to avoid this problem of non-uniqueness of 
binding energy curves, a random initial guess is used first 
for a configuration at small separation, in the knowledge 
that the result will be independent of the guess used. 
Then the centres computed at the previous, smaller sep- 
aration are used as the initial guess for the calculation 
at a larger separation. In this manner, a unique con- 
tinuous curve is obtained (labelled 'continuous' in Fig. [3] 
(bottom)). This is the approach that we adopt for all 
subsequent calculations in this paper. 

From the continuous curve, we obtain 3.97 A for the 
equilibrium separation and —0.28 kcal/mol for the bind- 
ing energy. This is in good agreement with the cou- 
pled cluster CCSD(T) calculations of Ref. [35|, which 
give 3.78 A and —0.28 kcal/mol, respectively, whereas 
revPBE without dispersion corrections gives 4.62 A and 
-0.04 kcal/mol. 
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Figure 3. Binding energy versus interatomic separation for 
the argon dimer, for varying relative orientations of the 
MLWF centres surrounding each atom (see Fig. [2j. Top 
panel: results obtained using the fragment method, in which 
the MLWF centres are calculated for a lone Ar atom and 
then translated and rotated to the second Ar atom. Bottom 
panel: results obtained using the true MLWF centres with 
various initial guesses for their positions. The curve labelled 
'continuous' is obtained by using the MLWF centres from a 
configuration at small separation as the initial guess for the 
centres at larger separations. In this way, the discontinuities 
in the curve are avoided and a unique curve is obtained (see 
text for details). 



C. Methane 

The methane dimer is a straightforward application of 
the Silvestrelli method: the positions of the MLWF cen- 
tres, which lie on the four tetrahedral C-H bonds of each 
CH4 molecule (see Fig. HJ, obey the same symmetries as 
the atomic positions, so there exists no arbitrariness of 
orientation. 

In Fig. [5j we compare to the results of both 
DFT+D and CCSD(T) calculations. Our geometries and 
CCSD(T) results were drawn from the Benchmark En- 
ergy and Geometry Database (BEGDB)^. 

The accuracy of Silvestrelli's method in the case of the 
methane dimer is good compared to CCSD(T): the for- 
mer gives an equilibrium separation of 3.66 A and bind- 




Figure 4. Illustration of the methane dimer. Carbon atoms 
are shown by large grey spheres, hydrogen by small white 
spheres, and the valence MLWF centres are shown by small 
pink spheres. 
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<y^> Semi-empirical DFT+D 
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Figure 5. Binding energy curves for the methane dimer with 
various methods. 



ing energy of —0.69 kcal/mol, and the latter 3.72 A and 
—0.53 kcal/mol, respectively. DFT+D is in somewhat 
worse agreement with CCSD(T), yielding 3.54 A and 
—0.76 kcal/mol, respectively. 



D. Ethene 

We now turn our attention to the ethene dimer, which 
includes a C-C double bond. Again we will compare 
results for the original and modified methods against 
CCSD(T) and DFT+D results. We have again used the 
geometries for each molecule taken from the BEGDB. 

To use Silvestrelli's original method in this case, we in- 
clude only the valence manifold in the creation of the ML- 
WFs, giving six MLWFs per molecule arranged as shown 
in Fig. [6] (left). In our modified method we use seven ML- 
WFs per molecule, with p-like, partly occupied orbitals 
on each carbon atom (Fig. [6] (right)). 

As seen in Fig. [7J neither the original Silvestrelli 
method (blue squares) nor DFT+D (red diamonds) re- 
produce the CCSD(T) values very accurately. By ex- 
panding the manifold of eigenstates used in the construc- 
tion of the MLWFs and applying our modified method 
to include partial MLWF occupancies and splitting of 
the p-like functions (see Sec. we find an excellent 
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Figure 6. Colours as in Fig. HJ Left: Ethene dimer with 
six MLWFs per molecule. Right: Ethene dimer with seven 
MLWFs per molecule. The centres of the p-like MLWFs are 
placed on the carbon atoms, but here we show the centres of 
the individual lobes of these p-like orbitals as calculated by 
our method. 
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Figure 7. Binding energy for an ethene dimer with various 
methods. 



agreement (black circles) with the CCSD(T) equilibrium 
values of 3.72 A for the separation and —1.51 kcal/mol 
for the binding energy; our method gives 3.73 A and 
— 1.60 kcal/mol, respectively; Silvestrelli's method gives 
3.83 A and -1.69 kcal/mol; DFT+D yields 3.55 A and 
-2.04 kcal/mol. 



E. Benzene 

For benzene, the valence states can be represented by 
15 doubly-occupied Wannier functions. The MLWF opti- 
mization procedure in this case therefore breaks the D&h 
symmetry of the benzene ring: the end result is that 



Figure 8. The three configurations used for the benzene dimer 
calculations: S (vertical displacement), PD (vertical and lat- 
eral displacement) and T (vertical displacement plus rotation 
in plane of one molecule), and the valence MLWF centres in 
each case (depicted by pink spheres). 



there are three C-C 'double' bonds and three C-C 'sin- 
gle' bonds in the MLWF representation. Those alternat- 
ing double and single C-C bonds represent a delocalised 
7r-bond around the ring. The double bonds are repre- 
sented by two centres located above and below the plane 
of the molecule, while the single bonds are represented by 
one centre on the bond. When two molecules are put in 
proximity (see Fig. [8j) and the vdW energy is calculated 
by Silvestrelli's method, the breaking of the symmetry 
affects the vdW energy in an arbitrary manner, depen- 
dent on how the two rings are aligned (i.e. whether the 
pairs of double bonds in adjacent molecules are aligned 
or anti-aligned). This alignment is defined by where the 
initial guesses for the centres of the Wannier functions 
are placed. 

The case of the benzene dimer therefore illustrates 
again the need to include the unoccupied antibonding 
states in the construction of the MLWFs: doing so in- 
creases the number of MLWFs to 18 and introduces par- 
tial occupancies, but restores the D&h symmetry of the 
system and also localises the MLWFs more. This then 
makes the vdW contribution independent of the initial 
guess for the Wannier function centres. 

We applied our implementation of the original Sil- 
vestrelli's method (with 15 MLWFs), and then our mod- 
ified method (with 18 MLWFs, partial occupancies and 
splitting of p-like states) to determine the binding energy 
as a function of displacement for three types of displace- 
ment (labelled S, PD, and T, illustrated in Fig. [8] of one 
of the molecules in the benzene dimer. We compare this 
to DFT+D and to the CCSD(T) calculations of Ref. Sll 
We note that we used the same bond lengths for C-C and 
C-H as Ref. l37l to within two decimal places, to construct 
perfectly symmetric benzene rings for our calculations. 

The binding energy curves for the various methods 
for the three configurations are shown in Fig. [9l Sil- 
vestrelli's method (blue squares) does not agree very 
well with CCSD(T) calculations, overestimating equilib- 
rium distances by 0.07-0.25 A (Table U) and overestimat- 
ing binding energies by 0.28-1.25 kcal/mol (Table HI). 
In particular, the dispersion curve obtained from Sil- 
vestrelli's method does not agree asymptotically with the 
DFT+D curve (red diamonds). In the T configuration 
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DFT-PBE 
Semi-empirical DFT+D 
b-b Silvestrelli (15 MLWFs) 
o-o This work (18 MLWFs) 
18 MLWFs (no p-splitting) 
* CCSD(T) (Janowski et al) 




DFT-PBE 
<^<> Semi-empirical DFT+D 
b-b Silvestrelli (15 MLWFs) 
G-O This work (18 MLWFs) 
* CCSD(T) (Janowski et al) 




DFT-PBE 

<y^> Semi-empirical DFT+D 
b-b Silvestrelli (15 MLWFs) 
G-O This work (18 MLWFs) 
* CCSD(T) (Janowski et al) 



Interdimer distance (A) 

Figure 9. Binding energy (kcal/mol) curves for the various 
methods for the benzene dimer in the S, PD and T config- 
urations (top, middle and bottom respectively). For the S 
configuration we also show the curve using 18 MLWFs per 
molecule if no p-splitting is used; in this case the method 
overbinds. CCSD(T) benchmark values are from Janowski et 



Silvestrelli's method performs better in terms of equilib- 
rium distance, binding energy and asymptotics as it can 
be seen in Fig. [9] (bottom). 

For the S configuration we also show the binding curve 
obtained if the anti-bonding states are included in the 
construction of the MLWFs, but splitting of the p-like 
states is not used (orange crosses); it is clear that in this 
case the method does not perform well, as replacing a p- 



Method 


S 


PD 


T 


Silvestrelli (15 MLWFs) 


4.01 


3.78 


5.06 


This work (18 MLWFs) 


3.89 


3.55 


4.88 


Semi-empirical DFT+D 


3.93 


3.58 


4.89 


CCSD(T) (Janowski et a&>) 


3.92 


3.53 


4.99 



Table I. Equilibrium distances in A for the benzene dimers in 
the three configurations (Fig. [8]) using the various methods. 
For all DFT calculations the PBE functional was used. 



Method 


S 


PD 


T 


Silvestrelli (15 MLWFs) 


-2.85 


-3.23 


-2.85 


This work (18 MLWFs) 


-1.47 


-2.31 


-2.64 


Semi-empirical DFT+D 


-1.38 


-2.11 


-2.87 


CCSD(T) (Janowski et a^I) 


-1.60 


-2.55 


-2.57 



Table II. Binding energies (kcal/mol) at equilibrium geometry 
for the benzene dimers in the three configurations (Fig. [8} 
using the various methods. For all DFT calculations the PBE 
functional was used. 



like orbital by an s-orbital is a very poor approximation. 

Our full modified method, including both the larger 
manifold and the splitting of p-like states (black circles 
in Fig. [9]), on the other hand, has excellent agreement 
in terms of equilibrium distances and binding energies 
with the DFT+D curves and the CCSD(T) values, for all 
three configurations, to within 0.05 A and 0.33 kcal/mol 
(Table HI and Hi]) ; the asymptotic behaviour of the energy 
is also better captured. 



F. H 2 Pc and CuPc 

To examine the difficulties encountered applying these 
methods to larger systems, we have investigated the ph- 
thalocyanine (H^Pc) dimer in the simplest configura- 
tion (S vertically displaced) first by applying Silvestrelli's 
method and then by applying our modifications it, and 
comparing the binding energy curve to one obtained us- 
ing DFT+D. The optimised MLWF centres for a single 
H2PC are shown in Fig. [10] (top). We see that as with the 
benzene molecule, there are alternating single and double 
MWLF centres on the C-C bonds of the six-membered 
rings, representing delocalised 7r-bonds. We also find, 
however, that using only the 93 valence MLWFs (186 va- 
lence electrons) is problematic, as a good representation 
of the electronic density of the system cannot be obtained 
in this way since this breaks the symmetry of the system, 
but most importantly it yields one lone MLWF of unre- 
alistically large spread (~2.5 A) located some distance 
from any atoms (Fig. [10] (top)). This is due to the fact 
that an odd number (93 MLWFs) is incompatible with 
the D2h symmetry of the molecule. 

Using a larger and even number of MLWFs (112 per 
molecule) we can restore this D2h symmetry of the 
molecule (Fig. [10] (bottom)) and represent the electronic 
density of the system in a way more compatible with its 
chemistry. When anti-bonding states are included, it is 



important to make a chemically intuitive initial guess for 
the centres and forms of the MLWFs. We make initial 
guesses as follows: we place p-like orbitals on the carbon 
atoms and s-like orbitals on every bond and p-like or- 
bitals on the hydrogenated nitrogens as well as two s-like 
orbitals on every non-hydrogenated nitrogen atom. In 
this way, we have partly occupied MLWFs that represent 
the 372 valence electrons of the dimer. The binding en- 
ergy curves obtained by using this representation and our 
modifications to Silvestrelli's method are shown in Fig.fTTl 
and compared to DFT+D. The binding energy obtained 
from our method is —23.63 kcal/mol and the equilibrium 
distance 3.58 A; with DFT+D we obtain -18.91 kcal/mol 
and 3.68 A. As for benzene, we see very good agreement 
with DFT+D; these values roughly agree with the stack- 
ing distance of crystalline H2PC (around 3.2-3.4 A) 38 . 
Silvestrelli's original method severely overbinds the dimer 
(giving a binding energy of —41 kcal/mol) because of the 
unphysically large spread of the lone MLWF that appears 
in the valence representation. This is due to the strong 
dependence of the vdW energy on the spreads (Eq. ([TB]) ). 

In the case of CuPc dimer (vertically displaced S con- 
figuration) we again do not use the valence manifold of 
390 MLWFs per dimer (195 MLWFs per molecule: 98 spin 
up and 97 spin down), but instead use a larger manifold 
of MLWFs. We note that the dimer configuration used 
here does not correspond to any phases CuPc is observed 
in experimentally, but was used for illustrative purposes 
as it is the simplest one. This is a spin-polarized system, 
so a different set of MLWFs is required for spin up / down 
electrons, yielding a total of 234 singly occupied MLWFs 
per molecule (117 for every spin channel). There are 10 
d-like MLWFs (five for every spin channel) centred on 
each copper atom, and s-like MLWFs on bonds and ni- 
trogens. The MLWFs corresponding to spin up and spin 
down electrons have essentially the same centres for the 
same bonds or atoms (Fig. [T2]) . 

In such cases, where some Wannier functions centres 
are very closely centred, it would be incorrect to consider 
them as separate fragments since this would violate the 
fundamental assumption of the ALL method, that it is 
valid for non-overlapping fragments only. This can be 
understood from the fact that Eq. ([T3]) is strongly non- 
linear, so adding the contributions of overlapping density 
fragments does not give the same result as summing the 
densities beforehand. As a result, Silvestrelli's method 
severely overbinds the dimer (~ —108 kcal/mol), demon- 
strating that the method breaks down for overlapping 
fragments. 

We alleviate this problem by amalgamating all the cen- 
tres and spreads of the closely placed MLWFs (in this 
case the d-like MLWFs on Cu) into one MLWF with a 
centre and spread given by the arithmetic mean of the 
closely placed MLWFs, and occupancies given by the sum 
of the separate MLWFs. The criterion for amalgamating 
MLWFs can be automated such that MLWFs less than 
a particular threshold distance apart are combined. In 
our case, we used a value of 0.1 A for this threshold, 
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Figure 10. Left: Phthalocyanine (H2PC) molecule and its va- 
lence MLWF centres. Hydrogen atoms are by small white 
spheres, carbon atoms by large grey spheres and nitrogen 
atoms by large blue spheres. The MLWF centres are shown 
by the small pink spheres. Using only the valence MLWFs 
does not give a satisfactory description of the system since 
it yields a lone MLWF of unphysically large spread (shown 
by large yellow sphere and labelled by the letter L). Right: 
H2PC molecule and its 112 MLWF centres, now including anti- 
bonding states. With this representation all the D 2 h symme- 
try of the ring is restored and a better chemical picture is 
given. There are s-like orbitals on every bond and the non- 
hydrogenated nitrogens, and p-like partly occupied orbitals 
on every carbon the two hydrogenated nitrogens (not shown 
here as these are located inside the corresponding atoms). 




Interdimer distance (A) 



Figure 11. Binding energy curves for H2PC dimer in the S 
configuration (vertically displaced) versus intermolecular dis- 
tance obtained with the various methods. 
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Figure 12. Copper phthalocyanine (CuPc) molecule and 
its 234 MLWF centres, again including anti-bonding states. 
Colours as in Fig. [TOj with copper shown by the large brown 
sphere in the centre. There are s- symmetry MLWFs on every 
bond and atom except for copper, p-like MLWFs on the car- 
bons and 5 d-symmetry MLWFs on the copper atom. Now 
there are no p-like orbitals on any nitrogen atom as for H2PC. 



System 


Cq (EhCLo) 




Silvestrelli 


This work 


MP2+AvdW 


pseudo-DOSD 


Argon 


92.4 


92.4 


76.1 


64.3 


Methane 


99.1 


99.1 


119 


130 


Ethene 


275 


261 


328 


300 


Benzene S 


2727 


1288 


2364 


1723 


Benzene PD 


2727 


1284 


2364 


1723 


Benzene T 


2769 


1262 


2364 


1723 



Table III. Effective intermolecular Cq coefficients. Dispersion- 
corrected MP2 (MP2+AvdW) and reference values are drawn 
from Ref. Eol . For the argon and methane dimers, our ap- 
proach is identical to the original method of Silvestrelli. The 
differences between the values reported in the first column 
(Silvestrelli) and those in Ref. are attributable to the differ- 
ent calculational details such as choice of exchange and corre- 
lation functional, simulation cell size and plane-wave energy 
cutoff. 



G. Intermolecular Cq coefficients 



1 1 1 1 1 1 



1 1 1 



3.6 3.8 4.0 4.2 4.4 4.6 4.8- 



DFT-PBE 
<y^> Semi-empirical DFT+D 
G-o This work (468 MLWFs) 




It is expedient to define effective intermolecular Cq co- 
efficients, 



Interdimer distance (A) 



Figure 13. Binding energy curves for the CuPc dimer in the S 
configuration (vertically displaced) obtained using the various 
methods. 



which had the desired effect of including the d-\ike or- 
bitals on Cu in the amalgamation procedure, while leav- 
ing all other MLWFs in the system unaffected. 

In Fig. [13] we compare the binding energy curves ob- 
tained using DFT+D to our modified method (now in- 
cluding the amalgamation of closely-overlapping ML- 
WFs) using a larger manifold of 468 MLWFs per dimer. 
This gives much more sensible results, with a binding en- 
ergy of —27.22 kcal/mol and an equilibrium separation 
of 3.57 A, in fair agreement with DFT+D, which gives 
—22.21 kcal/mol and 3.63 A, respectively. These val- 
ues are in reasonable agreement with those for H2PC (as 
obtained using our method above), and also with those 
obtained with other methods for other metal phthalo- 
cyanines (NiPc and MgPc calculated with the TS-vdW 
scheme in Ref. [§9| using the PBE funtional). 



^6eff — ~ ^ CQ n i, 



(24) 



where only intermolecular terms are summed over, i.e., n 
and / correspond to MLWFs on different molecules, and 
the factor of 1/2 accounts for double-counting. In Ta- 
ble IIII| we compare our values to those of the original 
method of Silvestrelli, benchmark dispersion-corrected 
MP2 calculations (MP2+AvdW) and reference results 
obtained using the Dipole Oscillator Strength Distribu- 
tion (DOSD) approach, given in the database of Ref. Ho|. 

As previously discussed in Ref. [2|, comparison with ref- 
erence values is made somewhat difficult by the fact that 
they are obtained by fitting to experimental data and 
hence also include higher-order terms (Cg, C10) in an 
effective manner. 

Taking the reference values as a benchmark, it can be 
seen from Table ITTT1 that, for the systems under consider- 
ation, there is no clear or systematic improvement in cal- 
culated effective Cq coefficients with our modifications to 
Silvestrelli's approach as compared to Silvestrelli's origi- 
nal approach: in the case of ethene the original method 
compares more favourably, while in the case of the ben- 
zene dimers our approach performs much better. In spite 
of this, however, it is worth noting that our approach (as 
shown earlier) significantly improves the values obtained 
for equilibrium separations and binding energies, as com- 
pared to CCSD(T), for all systems considered for which 
we have access to CCSD(T) results. 



H. Sensitivity to cutoff radius r c 

The sensitivity of the binding energy on the cutoff ra- 
dius r c in Eq. [Tfjlwas tested on the S configuration of the 
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3.40 3.60 3.80 4.00 4.20 4.40 4.60 4.80 5.00 

Interdimer distance (A) 

Figure 14. Binding energy curve for the benzene dimer in 
the S configuration for various values of r c using our modified 
method with 18 MLWFs per molecule. 

benzene dimer with 18 MLWFs per molecule (Fig. [T4|) . 
Even small changes of 1% in the cutoff radius result in 
significant changes in the binding energy curves, with 
the binding energy and equilibrium distance varying by 
6-8% and 0.2-0.8% respectively. For larger changes in r c , 
the method breaks down, as the energy changes are un- 
physically large. Although the cutoff radius is physically 
justified 7 , this strong dependence of the vdW correction 
on it is a weakness of the method. 



I. Approximations to the density 

In the original method of Silvestrelli, the KS density is 
approximated by replacing all real MLWFs with hydro- 
genic s wavefunctions (r) given by Eq. ([T5]) ; for the 
purpose of calulating the Cq coefficients, the electronic 
charge density of the system is, therefore, effectively ap- 
proximated as 

N occ 

p s {r)=Y;\<{*)?- (25) 

In the modified method presented here, in which the 
MLWFs are constructed using a manifold of the KS states 
beyond just the occupied orbitals, there are two levels 
of approximation to the charge density. First, the off- 
diagonal component Pod(?) is neglected from Eq. ([TO]) 
and, second, the "hydrogenic" approximation of the orig- 
inal approach is applied, whereby the disentangled Wan- 
nier functions, i^ ls (r), are replaced by hydrogenic or- 
bitals, w^(r), of the same center and spread. In our 
method, therefore, the density is approximated as 

N 

p dis (r) = J2J n>nto| 2 , (26) 

n=l 

where N is now the of total number of fragments, af- 
ter the splitting of p-like orbitals or amalgamation of co- 



| -2.40 ■ » I3» 
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Figure 15. Density profile on a plane parallel to a C-C bond 
(xz-plane) in a benzene molecule. Left: the original Kohn- 
Sham density p(r) from the plane-wave DFT calculation. 
Right: The off-diagonal component por>(r) of the density (see 
Eq. (|10p) when a disentangled manifold is used to construct 
A^dis = 18 MLWFs. Note the much-reduced scale compared 
to that of the total density. The units are A -3 . 




Figure 16. Density difference isosurface plots showing the dif- 
ference p(r) — p s (r) between the KS density and the approxi- 
mate "hydrogenic" density of the original Silvestrelli approach 
(Eq. J5SJ)). Left: cross-section through a "single" bond. Right: 
cross-section through a "double" bond. 

centric MLWFs has been performed. We consider each 
of these approximations in turn for a typical system, the 
benzene molecule. 

The XCrySDen 41 package was used to generate the 
isosurface plots referred to in this Section. 

In Fig. [15] we show density isosurface plots for the KS 
density p(r) (left) and the off-diagonal density Pod(t) 
(right), which emphasises that the latter is uniformly 
small in magnitude, comprising only a small fraction of 
the total density (~ 5 — 7%), as a result of the exponential 
localisation of the MLWFs. 

In Fig. [16] we show the difference between the KS den- 
sity p(r) and the "hydrogenic" approximation p s ( r ) of the 
original Silvestrelli method (Eq. ([25]) ) for two of the C- 
C bonds in benzene: on the left a "single" bond; on the 
right a "double" bond. These two bonds only differ be- 
cause of the symmetry-breaking inherent in the MLWF 
construction when just the valence states are used. We 
see that the density associated with the 7r-bond is not 
well represented in either case. 

Finally, in Fig. [T71 we show the difference between the 
KS density p(r) and that of our modified method Pdis(r) 
(Eq. ([26]) ) with 18 MLWFs obtained by disentanglement 
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Figure 17. Density difference isosurface plots, on the same 
plane as in Fig. [15] showing the difference p(r) — pdis(r) be- 
tween the KS density and the "hydrogenic" density of our 
method when a disentangled set of MLWFs is used (Eq. ([26]) ). 
Left: without splitting of p-like states; right: with splitting of 
p-like states into two s-like states. The mean difference with 
the KS density compared to the original Silvestrelli's method 
is reduced overall for both cases, but even more in the case of 
p-splitting. 



from a larger manifold. The left-hand plot is without 
splitting the p-like states, and the right-hand plot is with 
(as described in Sec. HIT]) . 

We see that while this introduces small regions where 
the density differs significantly (right at the MLWF cen- 
ters), everywhere else it is overall an improvement, pro- 
ducing a better representation of the density compared 
to the original Silvestrelli's method, especially in the case 
of p-splitting. 

In summary, discarding the off-diagonal component of 
the density (in the case of disentangled MLWFs) is a 
relatively minor approximation, and has a considerably 
smaller effect than approximating the density in various 
ways using hydrogenic orbit als, the latter being inher- 
ent to both our approach and the original approach of 
Silvestrelli. The maximum difference between the KS 
density p(r) and the density in our method is reduced by 
~ 23% and the minimum difference by ~ 5%, compared 
to the difference between the KS density and the density 
in Silvestrelli's method. 



J. Sensitivity to energy window E w - in 

To use our modifications to Silvestrelli's method, the 
disentanglement procedure has to be used in the con- 
struction of Wannier functions, as outlined in the Meth- 
ods section. Because including ever more high-energy 
plane- wave states inevitably allows extra variational free- 
dom in the construction of the MLWFs, we find that 
the precise values of the MLWF spreads are sensitive to 
the outer energy window used for the disentanglement. 
Specifically as E w [ n is increased, the MLWFs become 
more localised (their spreads decrease). As a result, the 
vdW energy is also affected by the choice of E w { n . 

In this work we have chosen throughout to estimate an 
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Figure 18. Binding energy curve for the benzene dimer in 
the S configuration for various values of a using our modified 
method with 18 MLWFs per molecule. 

appropriate energy window using 

^win = CLUMO + ^(eHOMO - eo) (27) 

where a is a factor used to scale the valence energy band- 
width. This is motivated by the idea that to enable us to 
restore the symmetry, we need to include the antibonding 
counterparts to the valence states, without including too 
large a number of irrelevant higher- lying unbound states. 
Eq. ([27|) is an attempt to estimate the range of energies 
spanned by these antibonding states. In Fig. [18] we show 
the dependence of the vdW binding energy curves for the 
benzene dimer in the S configuration on a. While there 
is considerable variation for too-small a, we find that for 
values beyond 0.4, the curves vary only a rather small 
amount with a; As long as a value of a around this value 
is chosen, it should yield reasonable results, suggesting 
the extra degree of empiricism introduced by this proce- 
dure is relatively limited in scale. The value of a was set 
to 0.4 in all the other calculations in this work. 

V. CONCLUSION 

We conclude that Silvestrelli's method is computation- 
ally efficient and very easy to implement for small sys- 
tems where initial guesses for the Wannier centres can be 
specified. However, there is a very strong dependence of 
the calculated vdW energy on the position and spread of 
the Wannier centres, and these are not always as unique 
as one might hope. Symmetry-breaking, often induced by 
considering only the valence manifold in the construction 
of the MLWFs, may introduce arbitary dependence on 
initial guesses in a way that significantly affect the vdW 
energy. We have shown that arbitrarily-broken symme- 
tries may often be restored by increasing the number 
of Wannier functions used and generating them with a 
suitably-chosen range of the conduction states as well as 
the valence states. This necessitates the inclusion of oc- 
cupancies in the formalism. We note that in cases where 
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no symmetries are restored when we use more MLWFs, 
as in the example of ethene, it is the better localisation of 
the MLWFs that may be responsible for improved vdW 
energies, since the method is based on pairwise summa- 
tion of well-separated fragments. 

Particularly, in cases with a larger number of Wan- 
nier functions, we have shown that the approximation 
implicit in replacing the true Wannier functions with 
hydrogenic s-orbitals may not always yield an accurate 
representation of the electronic density, and have shown 
how in cases where there is p-like symmetry, it is bet- 
ter to substitute the p- symmetry functions with two s- 
like functions. By considering the problems associated 
with applying these adapted methods to larger systems 
such as H2PC and CuPc, we have demonstrated that the 
approach is not necessarily a good candidate for study- 
ing larger systems, where specifying initial guesses for 
a large number of non-trivial MLWFs may be difficult; 
chemical insight for the form of these higher-lying states 
has to be employed, but becomes more difficult for even 
larger systems. In the case of copper phthalocyanine, we 
showed that MLWFs that are centred effectively at the 
same point (such as the five d-like MLWFs on each Cu 
atom) cannot be treated as separate fragments of density; 
they should instead be amalgamated into one fragment of 
density of an averaged centre and spread and summed oc- 
cupancies. The reason for this is that the method is valid 
only in the limit of well- separated fragments. Finally, we 
have demonstrated that there is also a strong dependence 
of the vdW energy on the cutoff radius used in the inte- 
gral of Eq. ([T6]h and although the value used is justified 
on physical grounds, it nevertheless represents something 
of an adjustable parameter with considerable influence 
on the results obtained. Overall, we conclude that while 
Silvestrelli's method suffers from several drawbacks, it 
can be made rather accurate once modifications are ap- 
plied to it (albeit with the introduction of further em- 
pirical character); these improvements, and Silvestrelli's 
method in general, however, may be less suitable for more 
structurally complex, large-scale systems, for which al- 
ternative methods that are more fully ab initio may be 
desirable. 
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